home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlahqr.f < prev    next >
Text File  |  1996-07-19  |  14KB  |  411 lines
  1.       SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
  2.      $                   ILOZ, IHIZ, Z, LDZ, INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     October 31, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       LOGICAL            WANTT, WANTZ
  11.       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DLAHQR is an auxiliary routine called by DHSEQR to update the
  21. *  eigenvalues and Schur decomposition already computed by DHSEQR, by
  22. *  dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
  23. *
  24. *  Arguments
  25. *  =========
  26. *
  27. *  WANTT   (input) LOGICAL
  28. *          = .TRUE. : the full Schur form T is required;
  29. *          = .FALSE.: only eigenvalues are required.
  30. *
  31. *  WANTZ   (input) LOGICAL
  32. *          = .TRUE. : the matrix of Schur vectors Z is required;
  33. *          = .FALSE.: Schur vectors are not required.
  34. *
  35. *  N       (input) INTEGER
  36. *          The order of the matrix H.  N >= 0.
  37. *
  38. *  ILO     (input) INTEGER
  39. *  IHI     (input) INTEGER
  40. *          It is assumed that H is already upper quasi-triangular in
  41. *          rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
  42. *          ILO = 1). DLAHQR works primarily with the Hessenberg
  43. *          submatrix in rows and columns ILO to IHI, but applies
  44. *          transformations to all of H if WANTT is .TRUE..
  45. *          1 <= ILO <= max(1,IHI); IHI <= N.
  46. *
  47. *  H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
  48. *          On entry, the upper Hessenberg matrix H.
  49. *          On exit, if WANTT is .TRUE., H is upper quasi-triangular in
  50. *          rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
  51. *          standard form. If WANTT is .FALSE., the contents of H are
  52. *          unspecified on exit.
  53. *
  54. *  LDH     (input) INTEGER
  55. *          The leading dimension of the array H. LDH >= max(1,N).
  56. *
  57. *  WR      (output) DOUBLE PRECISION array, dimension (N)
  58. *  WI      (output) DOUBLE PRECISION array, dimension (N)
  59. *          The real and imaginary parts, respectively, of the computed
  60. *          eigenvalues ILO to IHI are stored in the corresponding
  61. *          elements of WR and WI. If two eigenvalues are computed as a
  62. *          complex conjugate pair, they are stored in consecutive
  63. *          elements of WR and WI, say the i-th and (i+1)th, with
  64. *          WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
  65. *          eigenvalues are stored in the same order as on the diagonal
  66. *          of the Schur form returned in H, with WR(i) = H(i,i), and, if
  67. *          H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
  68. *          WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
  69. *
  70. *  ILOZ    (input) INTEGER
  71. *  IHIZ    (input) INTEGER
  72. *          Specify the rows of Z to which transformations must be
  73. *          applied if WANTZ is .TRUE..
  74. *          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
  75. *
  76. *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
  77. *          If WANTZ is .TRUE., on entry Z must contain the current
  78. *          matrix Z of transformations accumulated by DHSEQR, and on
  79. *          exit Z has been updated; transformations are applied only to
  80. *          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
  81. *          If WANTZ is .FALSE., Z is not referenced.
  82. *
  83. *  LDZ     (input) INTEGER
  84. *          The leading dimension of the array Z. LDZ >= max(1,N).
  85. *
  86. *  INFO    (output) INTEGER
  87. *          = 0: successful exit
  88. *          > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI
  89. *               in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
  90. *               elements i+1:ihi of WR and WI contain those eigenvalues
  91. *               which have been successfully computed.
  92. *
  93. *  =====================================================================
  94. *
  95. *     .. Parameters ..
  96.       DOUBLE PRECISION   ZERO, ONE
  97.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  98.       DOUBLE PRECISION   DAT1, DAT2
  99.       PARAMETER          ( DAT1 = 0.75D+0, DAT2 = -0.4375D+0 )
  100. *     ..
  101. *     .. Local Scalars ..
  102.       INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NR, NZ
  103.       DOUBLE PRECISION   CS, H00, H10, H11, H12, H21, H22, H33, H33S,
  104.      $                   H43H34, H44, H44S, OVFL, S, SMLNUM, SN, SUM,
  105.      $                   T1, T2, T3, TST1, ULP, UNFL, V1, V2, V3
  106. *     ..
  107. *     .. Local Arrays ..
  108.       DOUBLE PRECISION   V( 3 ), WORK( 1 )
  109. *     ..
  110. *     .. External Functions ..
  111.       DOUBLE PRECISION   DLAMCH, DLANHS
  112.       EXTERNAL           DLAMCH, DLANHS
  113. *     ..
  114. *     .. External Subroutines ..
  115.       EXTERNAL           DCOPY, DLABAD, DLANV2, DLARFG, DROT
  116. *     ..
  117. *     .. Intrinsic Functions ..
  118.       INTRINSIC          ABS, MAX, MIN
  119. *     ..
  120. *     .. Executable Statements ..
  121. *
  122.       INFO = 0
  123. *
  124. *     Quick return if possible
  125. *
  126.       IF( N.EQ.0 )
  127.      $   RETURN
  128.       IF( ILO.EQ.IHI ) THEN
  129.          WR( ILO ) = H( ILO, ILO )
  130.          WI( ILO ) = ZERO
  131.          RETURN
  132.       END IF
  133. *
  134.       NH = IHI - ILO + 1
  135.       NZ = IHIZ - ILOZ + 1
  136. *
  137. *     Set machine-dependent constants for the stopping criterion.
  138. *     If norm(H) <= sqrt(OVFL), overflow should not occur.
  139. *
  140.       UNFL = DLAMCH( 'Safe minimum' )
  141.       OVFL = ONE / UNFL
  142.       CALL DLABAD( UNFL, OVFL )
  143.       ULP = DLAMCH( 'Precision' )
  144.       SMLNUM = UNFL*( NH / ULP )
  145. *
  146. *     I1 and I2 are the indices of the first row and last column of H
  147. *     to which transformations must be applied. If eigenvalues only are
  148. *     being computed, I1 and I2 are set inside the main loop.
  149. *
  150.       IF( WANTT ) THEN
  151.          I1 = 1
  152.          I2 = N
  153.       END IF
  154. *
  155. *     ITN is the total number of QR iterations allowed.
  156. *
  157.       ITN = 30*NH
  158. *
  159. *     The main loop begins here. I is the loop index and decreases from
  160. *     IHI to ILO in steps of 1 or 2. Each iteration of the loop works
  161. *     with the active submatrix in rows and columns L to I.
  162. *     Eigenvalues I+1 to IHI have already converged. Either L = ILO or
  163. *     H(L,L-1) is negligible so that the matrix splits.
  164. *
  165.       I = IHI
  166.    10 CONTINUE
  167.       L = ILO
  168.       IF( I.LT.ILO )
  169.      $   GO TO 150
  170. *
  171. *     Perform QR iterations on rows and columns ILO to I until a
  172. *     submatrix of order 1 or 2 splits off at the bottom because a
  173. *     subdiagonal element has become negligible.
  174. *
  175.       DO 130 ITS = 0, ITN
  176. *
  177. *        Look for a single small subdiagonal element.
  178. *
  179.          DO 20 K = I, L + 1, -1
  180.             TST1 = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
  181.             IF( TST1.EQ.ZERO )
  182.      $         TST1 = DLANHS( '1', I-L+1, H( L, L ), LDH, WORK )
  183.             IF( ABS( H( K, K-1 ) ).LE.MAX( ULP*TST1, SMLNUM ) )
  184.      $         GO TO 30
  185.    20    CONTINUE
  186.    30    CONTINUE
  187.          L = K
  188.          IF( L.GT.ILO ) THEN
  189. *
  190. *           H(L,L-1) is negligible
  191. *
  192.             H( L, L-1 ) = ZERO
  193.          END IF
  194. *
  195. *        Exit from loop if a submatrix of order 1 or 2 has split off.
  196. *
  197.          IF( L.GE.I-1 )
  198.      $      GO TO 140
  199. *
  200. *        Now the active submatrix is in rows and columns L to I. If
  201. *        eigenvalues only are being computed, only the active submatrix
  202. *        need be transformed.
  203. *
  204.          IF( .NOT.WANTT ) THEN
  205.             I1 = L
  206.             I2 = I
  207.          END IF
  208. *
  209.          IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
  210. *
  211. *           Exceptional shift.
  212. *
  213.             S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
  214.             H44 = DAT1*S
  215.             H33 = H44
  216.             H43H34 = DAT2*S*S
  217.          ELSE
  218. *
  219. *           Prepare to use Wilkinson's double shift
  220. *
  221.             H44 = H( I, I )
  222.             H33 = H( I-1, I-1 )
  223.             H43H34 = H( I, I-1 )*H( I-1, I )
  224.          END IF
  225. *
  226. *        Look for two consecutive small subdiagonal elements.
  227. *
  228.          DO 40 M = I - 2, L, -1
  229. *
  230. *           Determine the effect of starting the double-shift QR
  231. *           iteration at row M, and see if this would make H(M,M-1)
  232. *           negligible.
  233. *
  234.             H11 = H( M, M )
  235.             H22 = H( M+1, M+1 )
  236.             H21 = H( M+1, M )
  237.             H12 = H( M, M+1 )
  238.             H44S = H44 - H11
  239.             H33S = H33 - H11
  240.             V1 = ( H33S*H44S-H43H34 ) / H21 + H12
  241.             V2 = H22 - H11 - H33S - H44S
  242.             V3 = H( M+2, M+1 )
  243.             S = ABS( V1 ) + ABS( V2 ) + ABS( V3 )
  244.             V1 = V1 / S
  245.             V2 = V2 / S
  246.             V3 = V3 / S
  247.             V( 1 ) = V1
  248.             V( 2 ) = V2
  249.             V( 3 ) = V3
  250.             IF( M.EQ.L )
  251.      $         GO TO 50
  252.             H00 = H( M-1, M-1 )
  253.             H10 = H( M, M-1 )
  254.             TST1 = ABS( V1 )*( ABS( H00 )+ABS( H11 )+ABS( H22 ) )
  255.             IF( ABS( H10 )*( ABS( V2 )+ABS( V3 ) ).LE.ULP*TST1 )
  256.      $         GO TO 50
  257.    40    CONTINUE
  258.    50    CONTINUE
  259. *
  260. *        Double-shift QR step
  261. *
  262.          DO 120 K = M, I - 1
  263. *
  264. *           The first iteration of this loop determines a reflection G
  265. *           from the vector V and applies it from left and right to H,
  266. *           thus creating a nonzero bulge below the subdiagonal.
  267. *
  268. *           Each subsequent iteration determines a reflection G to
  269. *           restore the Hessenberg form in the (K-1)th column, and thus
  270. *           chases the bulge one step toward the bottom of the active
  271. *           submatrix. NR is the order of G.
  272. *
  273.             NR = MIN( 3, I-K+1 )
  274.             IF( K.GT.M )
  275.      $         CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
  276.             CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
  277.             IF( K.GT.M ) THEN
  278.                H( K, K-1 ) = V( 1 )
  279.                H( K+1, K-1 ) = ZERO
  280.                IF( K.LT.I-1 )
  281.      $            H( K+2, K-1 ) = ZERO
  282.             ELSE IF( M.GT.L ) THEN
  283.                H( K, K-1 ) = -H( K, K-1 )
  284.             END IF
  285.             V2 = V( 2 )
  286.             T2 = T1*V2
  287.             IF( NR.EQ.3 ) THEN
  288.                V3 = V( 3 )
  289.                T3 = T1*V3
  290. *
  291. *              Apply G from the left to transform the rows of the matrix
  292. *              in columns K to I2.
  293. *
  294.                DO 60 J = K, I2
  295.                   SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
  296.                   H( K, J ) = H( K, J ) - SUM*T1
  297.                   H( K+1, J ) = H( K+1, J ) - SUM*T2
  298.                   H( K+2, J ) = H( K+2, J ) - SUM*T3
  299.    60          CONTINUE
  300. *
  301. *              Apply G from the right to transform the columns of the
  302. *              matrix in rows I1 to min(K+3,I).
  303. *
  304.                DO 70 J = I1, MIN( K+3, I )
  305.                   SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
  306.                   H( J, K ) = H( J, K ) - SUM*T1
  307.                   H( J, K+1 ) = H( J, K+1 ) - SUM*T2
  308.                   H( J, K+2 ) = H( J, K+2 ) - SUM*T3
  309.    70          CONTINUE
  310. *
  311.                IF( WANTZ ) THEN
  312. *
  313. *                 Accumulate transformations in the matrix Z
  314. *
  315.                   DO 80 J = ILOZ, IHIZ
  316.                      SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
  317.                      Z( J, K ) = Z( J, K ) - SUM*T1
  318.                      Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
  319.                      Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
  320.    80             CONTINUE
  321.                END IF
  322.             ELSE IF( NR.EQ.2 ) THEN
  323. *
  324. *              Apply G from the left to transform the rows of the matrix
  325. *              in columns K to I2.
  326. *
  327.                DO 90 J = K, I2
  328.                   SUM = H( K, J ) + V2*H( K+1, J )
  329.                   H( K, J ) = H( K, J ) - SUM*T1
  330.                   H( K+1, J ) = H( K+1, J ) - SUM*T2
  331.    90          CONTINUE
  332. *
  333. *              Apply G from the right to transform the columns of the
  334. *              matrix in rows I1 to min(K+3,I).
  335. *
  336.                DO 100 J = I1, I
  337.                   SUM = H( J, K ) + V2*H( J, K+1 )
  338.                   H( J, K ) = H( J, K ) - SUM*T1
  339.                   H( J, K+1 ) = H( J, K+1 ) - SUM*T2
  340.   100          CONTINUE
  341. *
  342.                IF( WANTZ ) THEN
  343. *
  344. *                 Accumulate transformations in the matrix Z
  345. *
  346.                   DO 110 J = ILOZ, IHIZ
  347.                      SUM = Z( J, K ) + V2*Z( J, K+1 )
  348.                      Z( J, K ) = Z( J, K ) - SUM*T1
  349.                      Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
  350.   110             CONTINUE
  351.                END IF
  352.             END IF
  353.   120    CONTINUE
  354. *
  355.   130 CONTINUE
  356. *
  357. *     Failure to converge in remaining number of iterations
  358. *
  359.       INFO = I
  360.       RETURN
  361. *
  362.   140 CONTINUE
  363. *
  364.       IF( L.EQ.I ) THEN
  365. *
  366. *        H(I,I-1) is negligible: one eigenvalue has converged.
  367. *
  368.          WR( I ) = H( I, I )
  369.          WI( I ) = ZERO
  370.       ELSE IF( L.EQ.I-1 ) THEN
  371. *
  372. *        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
  373. *
  374. *        Transform the 2-by-2 submatrix to standard Schur form,
  375. *        and compute and store the eigenvalues.
  376. *
  377.          CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
  378.      $                H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
  379.      $                CS, SN )
  380. *
  381.          IF( WANTT ) THEN
  382. *
  383. *           Apply the transformation to the rest of H.
  384. *
  385.             IF( I2.GT.I )
  386.      $         CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
  387.      $                    CS, SN )
  388.             CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
  389.          END IF
  390.          IF( WANTZ ) THEN
  391. *
  392. *           Apply the transformation to Z.
  393. *
  394.             CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
  395.          END IF
  396.       END IF
  397. *
  398. *     Decrement number of remaining iterations, and return to start of
  399. *     the main loop with new value of I.
  400. *
  401.       ITN = ITN - ITS
  402.       I = L - 1
  403.       GO TO 10
  404. *
  405.   150 CONTINUE
  406.       RETURN
  407. *
  408. *     End of DLAHQR
  409. *
  410.       END
  411.